

NASA Technical Memorandum 86830 


NAS A-TM-8 6830 19860005777 


Separated Transonic Airfoil Flow 
Calculations with a 
Nonequilibrium Turbulence Model 

L.S. King and D.A. Johnson 


, - ' „ i t 

November 1985 


■ i i 


imas/v 

National Aeronautics and 
Space Administration 



NF00072 



NASA Technical Memorandum 86830 


Separated Transonic Airfoil Flow 
Calculations with a 
Nonequilibrium Turbulence Model 

L. S. King, 

D. A. Johnson. Ames Research Center, Moffett Field, California 


November 1985 


NASA 

National Aeronautics and 
Space Administration 

Ames Research Center 

Moffett Field, California 94035 


/ ^ 

A lfo-/SZ¥7 



SUMMARY 


Navier-Stokes transonic airfoil calculations based on a recently developed 
nonequilibrium, turbulence closure model are presented for a supercritical airfoil 
section at transonic cruise conditions and for a conventional airfoil section at 
shock- induced stall conditions. Comparisons with experimental data are presented 
which show that this nonequilibrium closure model performs significantly better than 
the popular Baldwin-Lomax and Cebeci-Smith equilibrium algebraic models when there 
is boundary-layer separation that results from the inviscid-viscous interactions. 


INTRODUCTION 


The greatest limiting factor in the accurate numerical prediction of airfoil 
flows has been the lack of an adequate turbulence model. For airfoil flows with 
sufficiently mild pressure gradients and without separation, simple algebraic turbu- 
lence models (e.g., Cebeci and Smith, ref. 1) have been shown to yield good results 
(see e.g., refs. 2-4). As speeds are increased into the transonic regime and angles 
of attack are increased, however, adverse pressure gradients become stronger, and 
separated-flow regions make their appearance. For modern-day supercritical sec- 
tions, these problems are further complicated by the high surface curvatures on the 
rearward portion of the airfoil. Under these harsher conditions, the algebraic 
turbulence models do not work well. Characteristically, they overpredict pressure 
recovery on the rearward part of the airfoil. The separated-flow region is pre- 
dicted to be much thinner than that observed experimentally, and even the general 
shape of the velocity profile in the reverse-flow region does not agree with 
experiment. 

Recognizing the limitations in the simple algebraic models, considerable effort 
has been directed toward incorporating two-equation eddy-viscosity models into 
numerical prediction methods. Although these models are more complex, they are also 
more general than algebraic models. It was hoped that the greater generality would 
provide improved results for high-adverse-pressure-gradient and separated-flow 
situations. Indeed, improvements have been noted in predicting the experimentally 
observed rapid rise in skin friction downstream of reattachment, and mean-velocity 
profiles within the separated region are also in better qualitative agreement with 
experiment. Yet, for transonic cases with large separation, the prediction of shock 
position and surface pressures has been unsatisfactory (ref. 5). 

Recently, a turbulence closure model (ref. 6) designed specifically to treat 
two-dimensional, turbulent boundary layers with strong adverse pressure gradients 
and attendant separation has been developed. In this model, the influence of 
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history is modeled by using an ordinary differential equation. The equation, 
derived from the turbulence kinetic-energy equation, describes the streamwise 
development of the maximum turbulent shear stress. An eddy-viscosity distribution 
through the inner part of the layer is assumed which has as its velocity scale the 
maximum turbulent shear stress. In the outer part of the boundary layer, the eddy 
viscosity is treated as a free parameter which is adjusted to satisfy the Reynolds 
shear stress resulting from the ordinary differential equation. Because of this, 
the model is not simply an eddy-viscosity model, but contains features of a 
Reynolds-stress model. 

Results obtained with this new model incorporated into an inverse boundary- 
layer code were very encouraging (ref. 6). Subsequent to that work, further evalua- 
tions of this closure model were made (ref. 7) using a compressible Navier-Stokes 
method. The test flows in this latter study were those developed on the axisym- 
metric bump wind-tunnel model of reference 8. This wind-tunnel model was designed 
to produce inviscid-viscous interactions similar to those that develop on airfoils 
at transonic speeds. In the present work, the closure model has been incorporated 
into a Navier-Stokes airfoil code so that an assessment of its performance could be 
made on actual airfoil flows. Two airfoil test flows are considered. The first is 
the DSMA 671, supercritical airfoil section at test conditions selected to simulate 
transonic cruise. The second is the NACA 64A010 airfoil section at shock- induced 
stall conditions. The experimental data for these two test cases are reported in 
references 9 and 10, respectively. 


THE TURBULENCE MODEL 


The turbulence closure model as developed in reference 6 has as its basis three 
main observations: (1) that algebraic models such as the Cebeci-Smith model do 

quite well for attached flows in mild pressure gradients; (2) that they predict too 
rapid a rise in turbulent shear stress inside the boundary layer when high adverse- 
pressure gradients are encountered and, conversely, predict too rapid a decrease in 
this stress when the pressure gradients are relieved; and (3) that the inner mixing- 
length formulation of these models results m separated velocity profiles that are 
inconsistent with experiment. 

An algebraic eddy-viscosity distribution 


( 1 - exp - v fc /v. ) 

(1) 
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is used to describe the variation of the Reynolds shear stress normal to the shear 
layer under both attached and separated-flow conditions. In equation (2), D is a 
Van Dr iest -type near-wall damping term, y is the distance normal to the surface, 
and -u'v^ is the maximum Reynolds shear str ess divided by the local density at the 
given streamwise station. For convenience, -u'v^ will be referred to simply as the 

maximum Reynolds shear stress. In equation (3), a(x) is an unknown modeling param- 
eter; U e is the velocity at the edge of the shear layer; 6? is the incompressible 
boundary-layer displacement thickness; and 6 is the boundary-layer thickness. In 
the model, history effects (i.e., the slow response of the Reynolds shear stresses 
to local changes in the mean-velocity field) are taken into account through the 
parameter a(x). The streamwise distribution of this parameter is established in the 
solut ion procedure (details are given in ref. 7) to give a streamwise distribution 
of -u'v^ m the computed flow that satisfies an ordinary differential equation 

for -u'v^. This O.D.E., which is obtained from simplifications to the turbulence 
kinetic energy equation, can be found in references 6 and 7. 

The model described is designed for wall-bounded shear flows. In the wake, 
equation (3) was used to represent the eddy viscosities with separate displacement 
thicknesses used for the upper and lower parts of the wake. These displacement 
thicknesses were determined by integrating outward from the minimum wake-velocity 
locations. At the grid points corresponding to the location of minimum velocity, 
the eddy viscosity was taken to be the average o f th e eddy viscosities at the upper 
and lower adjoining grid points. In the wake, -u'v^ decreases with streamwise 

distance. As a result, a(x) will tend to increase from its value at the trailing 
edge. When a(x) was less than unity at the trailing edge (this generally will be 
the case), it was allowed to reach unity downstream but not allowed to exceed 
unity. If it were greater than 1 at the trailing edge (which would be more likely 
to occur on the lower airfoil surface), it was set equal to 1 in the wake. 

Although, initially the eddy viscosities in the upper and lower parts of the wake 
were unequal, farther downstream the wake was observed to develop symmetry causing 
the eddy-viscosity distributions also to become symmetric. 

In the very-far wake, equation (3) with o(x) = 1 underestimates the eddy 
viscosities by approximately a factor of 4 (ref. 11). However, it seems reasonable 
that the flow about the airfoil itself should not be sensitive to the very-far-wake 
development. The treatment of the wake is acknowledged to be very approximate. But 
m our opinion, the resultant airfoil surface pressures are much more sensitive to 
the turbulent Reynolds shear-stress development along the airfoil surface than to 
the rate of decay of these stresses in the wake. The favorable results obtained 
with this simple wake model support this opinion. 

In the two sets of calculations to be presented, a different velocity scale was 
used for the Van Driest damping term m equation (2). In the DSMA 671 calculations, 
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the conventional friction velocity u^ = (t /p ) was used. But in the NACA 
64A010 calculations, steady-state solutions could not be obtained with either the 
Cebeci-Smith model or the present model when u T was used. As in the supersonic 
compression corner study of reference 12, this problem was alleviated by using 
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(-u'v^) as the near-wall damping velocity scale. In both sets of calculations, 
the Van Driest constant, A + , was taken to be 15. 


THE NUMERICAL PROCEDURE 


The basic numerical method used in the present investigation is that due to 
Steger (ref. 13) for the Reynolds-averaged, time-dependent, compressible Navier- 
Stokes equations. In Steger 's method, the governing equations are transformed to a 
generalized body-fitted coordinate system and solved with the second-order-accurate, 
factorized implicit algorithm of Beam and Warming (refs. 14 and 15). In the origi- 
nal version of the code, viscous terms in the streamwise direction were neglected, 
resulting in the so-called "thin-layer approximation." The present version of the 
code contains all viscous terms, with an implementation similar to that of Degani 
(ref. 16). To account for wind-tunnel wall interference effects so that direct 
comparisons with wind-tunnel data may be made, the Steger code has been further 
modified by incorporating a pressure boundary condition along the upper and lower 
boundaries. Details of this modification are presented in earlier papers by King 
and Johnson (refs. 2 and 3). 

Mesh generation was accomplished using a Poisson solver similar to that of 
Thompson et al. (ref. 17) as modified by Steger and Sorenson (ref. 18). The mesh 
code produces a "wraparound," or C-mesh, and is coincident with user-prescribed 
points on the boundaries, that is, the airfoil surface and the outer computational 
boundary. With the Steger-Sorenson modification, orthogonality at the airfoil 
surface and concentration of coordinate lines near the surface may be controlled to 
resolve the turbulent boundary layer. In the DSMA 671 calculations, the mesh was 
composed of 139 points in the wraparound direction and 50 points m the direction 
away from the airfoil. For the NACA 64A010 calculations, the number of points in 
the wraparound direction was increased to 159 to ensure adequate resolution of the 
strong shock-wave/turbulent-boundary-layer interaction of this test case. With this 
more refined mesh, the streamwise spacing in the vicinity of the shock was approxi- 
mately 0.01 chord. The meshes were constructed with the first coordinate line off 
the airfoil at a normal distance of 2x10”^ chords from the surface. This distance 
corresponds roughly to a value of y + of 2, with approximately 20 points in the 
turbulent boundary layer near the airfoil midchord. 

The upper and lower mesh boundaries at which measured static pressures were 
applied as boundary conditions were located at ±1.125 and ±1.0 chord for the 
DSMA 671 and NACA 64A010 sections, respectively. The upstream boundary was located 
2 chords upstream of the airfoils. Since no data were available at the upstream 
boundary, an approximate set of upstream boundary conditions had to be created. For 
the DSMA 671 airfoil, these were obtained from inviscid free-air calculations (see 
ref. 2). In the NACA 64A010 calculations, uniform flow at the upstream boundary was 
assumed. 
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RESULTS AND DISCUSSION 


The results obtained for the DSMA 671 airfoil section at test conditions 
selected to simulate transonic cruise will first be presented. The test conditions 
were = 0.72, a = 4.32°, and Re c = 2.67*10 . Transition was initiated on the 
upper and lower surfaces at x/c = 0.17, which corresponds to the transition strip 
locations of the experiment. The aspect ratio of the airfoil was 3, and the tunnel 
half-height to chord was 1.5. 

The significance of wall effects at these test conditions is illustrated in 
figure 1 , which shows surface-pressure distributions obtained with the Baldwin-Lomax 
closure model (ref. 19), using free-air and pressure-boundary-conditions. The free- 
air solution was obtained with the upper and lower boundaries far removed from the 
airfoil. In order to make comparisons of computations and experiment, it is obvious 
that the pressure boundary condition (PBC) is necessary. This is also true for the 
NACA 64A010 test case. 

In figure 2, the surface pressures predicted with the present model and with 
the Baldwin-Lomax model, both using the PBC, are compared with the experimental 
results. As can be seen, the differences between the solutions of the present model 
and of the Baldwin-Lomax model are small on the forward portion of the airfoil. The 
small disagreement with experiment on the forward portion of the airfoil is believed 
to be a result of the upstream boundary conditions employed. It appears that these 
boundary conditions created an effective angle of attack that was slightly too 
large. On the rear portion of the airfoil, high adverse pressure gradients exist 
because of the high curvature. Separation is predicted to occur slightly forward of 
that observed in the experiment, at 0.95 chord instead of 0.98 chord. The inset of 
figure 2, showing the trailing-edge region in greater detail, illustrates that the 
solution based on the present model is m excellent agreement with the data in this 
high gradient and separated-flow region. On the other hand, the Baldwin-Lomax model 
overpredicts pressure recovery in the separated region. 

Predicted upper-surface boundary-layer velocity profiles are shown m figure 3 
along with the experimental data. Four stations are shown, from x/c = 0.63, down- 
stream of the shock, to x/c = 0.99, m the separated region. At x/c = 0.63, the 
pressure gradient is small, and the two results agree well with each other and with 
the data. Proceeding farther downstream, high adverse pressure gradients are 
encountered, with the result that increasingly larger differences between the two 
solutions can be seen. At x/c = 0.99, the flow is separated. At this streamwise 
station, the present model results are in good agreement with the experimental 
data. Such is not the case for the Baldwin-Lomax model results. With this closure 
model, the momentum loss incurred by the boundary layer at this station is substan- 
tially underpredicted, and the shape of the separated profile is not in agreement 
with the experimental results. 

In figure 4, the predicted and measured Reynolds shear-stress profiles are 
compared at the streamwise stations of the mean-velocity profiles of figure 3. As 
recommended m reference 7, these results are compared m shear-layer coordinates as 


5 



defined by the direction of the flow at the location of maximum Reynolds shear 
stress. The predicted stresses are lower than those measured, but there was concern 
m the experiment that the measured st resses may have been in error on t he high 
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side. This concern arose because v' was measured to be as large as u' across 
the boundary layer at x/c = 0.63 and 0.75. Skin-friction determinations from law- 
of-the-wall fits of the mean-velocity data at these two stations where the stream- 
wise pressure gradients were small agree better with the predicted shear stresses 
than with the measured shear stresses, which further suggests that the measured 
values of -u'v' may have been high. Notice from this figure the different rates 
at which the shear stresses are predicted to inc rease near the trailing edge for the 
two closure models. The slower growth m -u'v' predicted by the present model 
accounts for the larger momentum losses and lower pressure recoveries predicted by 
this model. 

As a result of improved velocity predictions with the present model, flow 
angles in the near-wake are predicted much better than with the Baldwin-Lomax 
model. The flow angle, 0 = arctan v/u, is shown m figure 5 along a vertical line 
in the near-wake 0.02 chord downstream of the trailing edge. In this figure, 
y = 0 corresponds to the vertical position of the airfoil trailing edge. Both the 
data and the calculations based on the present model show a large jump in flow angle 
where the boundary-layer flows from the upper and lower surfaces meet. 

Velocity profiles in the wake of the supercritical section are shown m 
figure 6. Because the velocity profiles of the present model were in substantial 
agreement with the data at the airfoil trailing edge, results in the near-wake also 
agree well with the data. Farther from the trailing edge, the results from the 
present model show a larger velocity defect than the data. This is as expected, 
since the simple wake model employed in the calculation results in eddy-viscosity 
values which are too low in the far-wake. In contrast to the Baldwin-Lomax model, 
however, the present model does predict wake position very well. 

Results are presented for the NACA 64A010 conventional airfoil m fig- 
ures 7-10. The test conditions were = 0.8, a = 6.2°, and Re c = 2x10 . Transi- 
tion was initiated along the upper and lower surfaces at x/c = 0.17, the transition 
strip locations of the experiment. The aspect ratio was 4 and the tunnel half- 
height to chord was 2. At these test conditions, the upper-surface boundary layer 
separates at the shock, and the flow remains separated into the wake. The thick, 
detached shear layer that develops downstream of the shock is clearly evident in the 
infinite-fringe mterferogram presented m figure 7(a). The fringes represent lines 
of constant density, which in the inviscid flow regions correspond very nearly to 
lines of constant Mach number. In figure 7(b), a Mach contour plot of the solution 
obtained with the present model is presented. 

In figure 8, surface-pressure predictions obtained with the present model and 
with the Cebeci-Smith model are compared with the experimental results. A steady- 
state solution could not be obtained for this test flow with the Baldwin-Lomax 
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model, so the Cebeci-Smith model (with (-u'v^) used in the near-wall damping 
expression) was run instead for comparison purposes. These two algebraic models 
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are, in theory, equivalent (ref. 19). As evident from figure 8, the present model 
does much better at predicting the resultant pressure distribution for this shock- 
induced stall condition. The Cebeci-Smith model significantly underpredicts the 
viscous displacement effects downstream of the shock, as seen from the mean-velocity 
profile comparisons presented in figure 9. The present model, on the other hand, 
does quite well at predicting the thick, detached shear layer that is generated by 
the shock wave farther upstream. Shown alongside the mean-velocity profiles are the 
predicted and measured Reynolds shear stresses. The calculations and the measure- 
ments are compared in shear-layer coordinates as defined by the direction of the 
flow at the location of maximum Reynolds shear stress. At the two upstream measure- 
ment stations, the shear stresses predicted by the present model are in quite good 
agreement with the experiment. The disagreement between the predicted and measured 
stresses just downstream of the trailing edge (x/c = 1.02) is the result of the 
approximate wake model. Near the minimum-velocity location, the eddy viscosities 
are too large, a result of equation (2) being abruptly dropped from the eddy- 
viscosity relationship immediately downstream of the trailing edge. This deficiency 
could easily be corrected, for example, by modeling the inner boundary-layer eddy- 
viscosities approach to the larger wake values to be a function of the wake-velocity 
deficit. 


Shown in figure 10 is the development of -u'v^ with streamwise distance as 

predicted with the present model and the Cebeci-Smith model. With the present 

model, a less rapid rise in -u'v^ is predicted to occur at the shock with -iu ' 

continuing to grow downstream as a detached shear layer is formed. With the Cebeci- 

Smith model, however, -u'v' attained its maximum value at the shock and then mono- 
’ ’ m 

tonically decreased downstream. The decay in -u'v^ downstream of the shock indi- 
cates that the computed boundary layer in this case had little difficulty negotiat- 
ing the larger adverse pressure gradients predicted on the aft section of the 
airfoil. 


CONCLUDING REMARKS 


A nonequilibrium turbulence model has been incorporated into a Navier-Stokes 
airfoil code. Solutions obtained for two different airfoil flows, one with a small 
trailing-edge separation bubble and the other with a large, shock-induced, detached 
shear layer, have shown that results with this nonequilibrium model are clearly 
superior to those obtained with equilibrium models like those of Baldwin and Lomax 
and Cebeci and Smith. Airfoil surface pressures and boundary-layer and wake- 
velocity profiles are in much better agreement with experimental data when the 
nonequilibrium model is used. 
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FREE-AIR, BALDWIN- LOMAX MODEL 

PBC, BALDWIN-LOMAX MODEL 

O EXPERIMENT 



Figure 1.- Comparison of PBC and free-air surface-pressure with experimental results 
for DSMA 671 supercritical airfoil section: M ro = 0.72, a = 4.32°, 

Re c = 2.67x10 . 
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CLOSURE MODEL 



Figure 2.- Comparison of computed and experimental surface pressures: 
supercritical airfoil, M m = 0.72, a = 4.32°, Re c = 2.67* 10^. 


DSMA 671 


O EXPERIMENT 

PRESENT MODEL 

BALDWIN - LOMAX MODEL 




Figure 3.- Upper-surface mean-velocity profiles for DSMA 671 section. 
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PRESENT MODEL 
BALDWIN-LOMAX MODEL 
EXPERIMENT 



Figure 4.- Upper-surface Reynolds shear-stress profiles for DSMA 671 section. 


PRESENT MODEL 

BALDWIN-LOMAX MODEL 

O EXPERIMENT 



0 = arctan v/u 


Figure 5.- Flow angles at x/c = 1.02 for DSMA 671 section. 
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O EXPERIMENT 

— PRESENT MODEL 

- - BALDWIN - LOMAX MODEL 



CONTOUR 
INCREMENT = 0.05 



(a) Infinite-fringe interferogram. (b) Mach contour plot of computation 

base on present model. 

Figure 7.- NACA 64A010 airfoil section at shock-induced stall conditions: 

M m = 0.8, a = 6.2°, Re c = 2x10 6 . 
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Figure 8.- Comparison of computed and experimental surface pressures for NACA 
64A010 airfoil at shock-induced stall conditions: M = 0.8, a = 6.2°, Re. = 2x10 

co 1 7 C 


□ EXPERIMENT PRESENT MODEL * CEBECI-SMITH MODEL 


x/c = 0.67 x/c = 0.83 



x/c = 1.02 



Figure 9.- Upper-surface and near-wake mean-velocity and Reynolds shear-stress 

profiles for NACA 64A010 section. 
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Figure 10.- Upper-surface maximum Reynolds shear-stress development with streamwise 

distance for NACA 64A010 section. 
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